SUBROUTINE G_RADFLX
!-----------------------------------------------------------------------
!
! Definition von G:  Strahlungsflussgleichung
!

      use primvar,  only : G, X, XA, MR, MUr, MD, MJ, MH, Gnorm
      use physco,   only : z1, z2, z3, z4, z12, z14, StBolz, pi, clight, offset
      use config,   only : np
      use global,   only : tst
      use geomvar,  only : V_flux, V_vol, V_volA, S_vol, V_vol
      use matvar,   only : OPAros
      use eddvar,   only : f_edd, mu_r
      use zvar,     only : rho0
      use advecvar, only : H_adv
      use RBvar,    only : L_star


      implicit none

      integer       :: i
      logical, save :: init = .true.

      if (init) then
         write(66,"(a)") "G_radflx.f90:   rechne voll implizit;   Innenrand H=0 ;  Aussenrand: H=0"
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      G(MH,1)     =  X(MH,1) - X(MH,1+1)
      Gnorm(MH,1) = max( abs(X(MH,1)), offset )
      G(MH,2)     =  X(MH,2) - X(MH,2+1)
      Gnorm(MH,2) = max( abs(X(MH,2)), offset )

   ! aeussere Pseudozellen: i=np & i=np-1
      G(MH,np)       =  X(MH,np) - X(MH,np-1)
      Gnorm(MH,np)   = max( abs(X(MH,np)), offset )
      G(MH,np-1)     =  X(MH,np-1) - X(MH,np-1-1)
      Gnorm(MH,np-1) = max( abs(X(MH,np-1)), offset )

   ! innere Randbedingung: i=3
      !~ G(MH,i)     =  X(MH,i) - L_star / ( X(MR,i) * z4*pi )**2 ! Einstrahlung
      G(MH,3)        =  X(MH,3)                            !    fixer Randwert - hier 0
      Gnorm(MH,3)    = max( abs(X(MH,3)), offset )
   ! aeussere Randbedingung: i=np-2
      !~ G(MH,i)     =     X(MH,i) - z12 * ( S_vol(i) * X(MJ,i) + S_vol(i-1) * X(MJ,i-1) ) / V_vol(i) * mu_r ! Abstrahlungsrandbedingung
      G(MH,np-2)     =  X(MH,np-2)                         !    fixer Randwert - hier 0
      Gnorm(MH,np-2) = max( abs(X(MH,np-2)), offset )


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=4,np-3

         G(MH,i) = V_vol(i) * X(MH,i) - V_volA(i) * XA(MH,i)                                                         & ! temporal diff.
                 + V_flux(i) * H_adv(i) - V_flux(i-1) * H_adv(i-1)                                                   & ! advection
                 + clight * z2*pi * X(MR,i) * (f_edd(i) * X(MJ,i) - f_edd(i-1) * X(MJ,i-1)) * tst                    & ! grad(K) term
                 + clight * z14 * ( X(MJ,i)   * ( z3*f_edd(i)   - z1 ) * S_vol(i)                                    & ! (3K-J) term
                                  + X(MJ,i-1) * ( z3*f_edd(i-1) - z1 ) * S_vol(i-1) ) / X(MR,i) * tst                & ! ...
                 + X(MH,i) * pi * X(MR,i) * ( X(MUr,i+1) - X(MUr,i-1) ) * tst                                        & ! H div(u) term
                 + clight * z12*( OPAros(i)*rho0(i)*S_vol(i) + OPAros(i-1)*rho0(i-1)*S_vol(i-1) ) * X(MH,i) * tst      ! kappa rho H

         Gnorm(MH,i) = max( &
              abs( V_vol(i) * X(MH,i) ),                                                                             & ! cell content
              abs( V_flux(i) * H_adv(i) ), abs( V_flux(i-1) * H_adv(i-1) ),                                          & ! advection
              abs( clight * z2*pi * X(MR,i) * (f_edd(i) * X(MJ,i) - f_edd(i-1) * X(MJ,i-1)) * tst ),                 & ! grad(K) term
              abs( clight * z14 * ( X(MJ,i)   * ( z3*f_edd(i)   - z1 ) * S_vol(i)                                    & ! (3K-J) term
                                  + X(MJ,i-1) * ( z3*f_edd(i-1) - z1 ) * S_vol(i-1) ) / X(MR,i) * tst ),             & ! ...
              abs( X(MH,i) * pi * X(MR,i) * ( X(MUr,i+1) - X(MUr,i-1) ) * tst ),                                     & ! H div(u) term
              abs( clight * z12*( OPAros(i)*X(MD,i)*S_vol(i) + OPAros(i-1)*X(MD,i-1)*S_vol(i-1) ) * X(MH,i) * tst )  & ! kappa rho H
                          )

      end do


END SUBROUTINE G_RADFLX



